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ABSTRACT. An algorithm is given to efficiently compute L-functions 
with large conductor in a restricted range of the critical strip. Ex- 
amples are included for about 21000 dihedral Galois representa- 
tions with conductor near 10 7 . The data shows good agreement 
with a symplectic random matrix model. 



1. Introduction. 

In |[12J, Odlyzko and Schonhage developed an algorithm to com- 
pute the Riemann zeta function ((s) efficiently for values of s very 
high up in the critical strip. Their method depends on precompu- 
tation of Taylor series expansions of ((s) at regularly spaced points, 
which in turn can be done efficiently by a clever application of the 
Fast Fourier Transform. Rumely later implemented a version of this 
for Dirichlet L-f unctions in il6l . 

In analytic number theory it is often the case that there is a sym- 
metry between what one can prove for large values of t = Im(s) in 
the critical strip, and what one can prove for L-functions with large 
conductor q. With this in mind, it seems reasonable to try to find 
an algorithm to efficiently compute values of L-functions with very 
large conductor. 

Roughly speaking, we want to view the L-function as a Mellin 
transform of an automorphic form /, and split the integral at the 
symmetry point. This gives the extended L-function as an infinite 
sum 

A(s, f) = J2 a{n){G(s, 2nn/q 1 ' 2 ) + G{1 - s, 2nn/q 1 / 2 )} 

n 

of values of an incomplete Gamma function 

G(s, x) = x~ s T(s, x) = [ exp(—yx)y s —. 

Ji y 
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It was already observed in ilZl that the rapid decay of G(s,x) as 
x — > oo means that we can truncate the infinite series at about 0(q 1 ^ 2 ) 
terms. On the other hand, the exponential decay as t = Im(s) in- 
creases will cause a loss of precision if we want to compute L-function 
values for large t. This can be fixed by moving the contour integral 
so that the x parameter has complex values, as observed in |10| and 
|8J and implemented in iTTSlT . 

The algorithm is based on computing Taylor series expansions of 
the function G(s, x) in the second variable. Differentiation under the 
integral sign and integration by parts gives a recursion which al- 
lows the computations of derivatives of G{s,x) at very little cost. 
The rapid decay of G(s, x) as x increases will imply that, instead of 
taking equally spaced points for the centers of the Taylor expansions, 
the sequence of points can grow exponentially. This, in turn, will im- 
ply that we need very few Taylor expansions to compute efficiently, 
in fact 0(log(g)) expansions with 0(log(g)) terms in each. 

However, too much change in t will require that we increase the 
phase in the x parameters, and then all Taylor expansions must be 
recomputed. For that reason, we will only consider \t\ < 1 in this 
paper, although the method will generalize to any bounded region 
in the critical strip. For further simplification we will consider only 
s = 1/2 + it on the critical line, although this, too, is not crucial. 

If we want to get the most benefit of the precomputations it makes 
sense to compute zeros of lots of L-functions with the same conduc- 
tor simultaneously. Here the power of the Fast Fourier Transform 
can again play a role, if all of the L-functions arise from the same 
abelian group structure. In this paper we will focus primarily on 
L-functions attached to characters <fi on the ideal class group of a 
complex quadratic field with discriminant —q. The corresponding 
automorphic forms are weight one forms, which are cusp forms if 
the character is not a genus character. The corresponding Galois 
representations are dihedral. (The final section of the paper, how- 
ever, includes an example of an elliptic curve L-function with large 
conductor, and of an L-function for a quadratic Dirichlet character 
for large discriminant.) 

With this in mind we now explain the algorithm in a little more 
detail. It will be convenient to use the correspondence between ideal 
classes and binary quadratic forms Q of discriminant —q. Corre- 
sponding to a character <p we have a theta function 

Q(z,(f>) = 4>{Q) r Q( n ) exp(27rmz) = f ' yV^,(n) exp(27imz). 

[Q] n>0 n>0 
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Then 

A(s,4>) d = (q 1/2 /2nyr(s)L(sA) = 

^r>(n) {G(s, 27rn/g 1 / 2 ) + G(l - s, 2 7 rn/ g 1 / 2 )} . 

n>0 

(See \5, §12.4] or (6l §22.3].) For purposes of locating zeros on the 
critical line it is better to compute the Hardy function Z{t, <fi) defined 
for t > by 

6(t, 0) = t log(g 1 / 2 /(27r)) + arg(r(l/2 + it)) 
and, with s — 1/2 + it, 

Z(t, 0) = exp(i9(t, <f>))L(s, 0) 

= exp^(t,0))(27r/g 1/2 )T(s)- 1 x 

r<(,{n) {G(s, 2<im/q 112 ) + G(l - s, 2nn/q 1 ' 2 )} . 

We will truncate the infinite series after N terms, and arrange the n < 
N into intervals Ij, j = 1, 2, ... T centered at points Xj and with width 
Aj. For n in the interval Ij we compute a Taylor series expansion, 
truncated to B terms 

B 

G(s, 27in/q 1/2 ) + G(l - s, 2nn/q 1/2 ) « ^ Gj, k (s)(27rn/q 1/2 - Xj ) k , 

k=0 

where 

(1) G j>k (s) = [G^ k \s,Xj) + C7 (fc) (l - s,Xj)} jk\. 
This gives Z(t, 0) as 

(2) Z{t,(j>) w exp^(t,0))(27r/g 1/2 )T(s)- 1 x 

j=l fc=0 ne/j ■? 

Here we have inserted canceling terms A*, Aj fc to control the small 
size of Gj t k(s) and the large size of (2Tcn/q 1 ^ 2 — a^-)*. 

The point, then, is that the inner sum over n E Ij is independent 
of s, and so can be done as a precomputation. Of course, this is 
done separately for each character 0, but the coefficients r^n) can 
be computed for all characters very efficiently via the Fast Fourier 
Transform once the representation numbers tq (n) are known. 

If we want to use this method for other kinds of L-functions, for 
example, the L-function of an elliptic curve, or a quadratic Dirichlet 
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character, we forgo the advantage of using the Fast Fourier Trans- 
form. However, we see below we can choose the parameters Xj to 
be an integer times (a constant analogous to) 2n/q 1 ^ 2 , so that the pre- 
computation can even be done with integer arithmetic when the co- 
efficients r(n) are integers. The A* 1 terms are simply omitted, and 
the (analog of the) (2n/q 1 /' 2 ) k term is factored out of the inner sum. 
The author would like to thank Kimberly Hopkins for help with this 
idea. 

The remaining sections of the paper are as follows: 

2. Review of properties of G(s, x). 

3. Truncation of the L-series after N terms. 

4. Arrangement of the Taylor expansions. 

The T intervals Ij = (xj — Aj/2, Xj + A^/2) are determined. 

5. Truncation of the Taylor expansions after B terms. 

6. Implementation and examples. 



2. Review of properties of G(s, x). 

As above we define for x > and s e C 

f°° dy 
G(s,x)=x s T(s,x) = / exp(— xy)y s — . 

Ji y 

Differentiating with respect to x under the integral we see that 

dy 

exp(—xy)y s+1 — = — G(s + 1, x). 

y 

Integration by parts, on the other hand, gives 

(4) G(s + l,x) = 6XP( ~ X) + -G(s, x). 

X X 

Equations © and © give a nice recursive relation for all the deriva- 
tives G^ k \s, x) in terms of G(s, x). 

We will need bounds for G(s, x) and its derivatives. For s in the 
critical strip, that is < Re(s) < 1, we have 

/oo 
exp(-xy)y Re ^- l dy. 

Change the variables by u = y — 1 to get 

POO 

\G(s,x)\ <exp(-x) / exp(-xu)(u + lf^^du. 
Jo 



(3) T GM=-I 
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Now u + 1 > 1 and by hypothesis, Re(s) - 1 < 0, so (u + l) Re («)-i < i, 
and therefore 

(5) |G(s, x)\ < exp(— x) [ exp(—xu)du = 



To estimate the derivatives G^ k \s, x) = (—l) k G(s + k, x) this method 
does not apply, since we only have Re(s) + k — 1 < k. Instead we can 
use the Cauchy formula for derivatives 

\GW(S,X)\ Mr 

k\ ~ R k ' 

where M R is a bound for G(s,w) with |w — x| = R. An estimate 
similar to © shows that |c7(s,w)| < exp(— Re(w))/Re(w), which a 
calculus argument shows is maximized at w = x — R, that is, 

|C?( fc )( g ,x)| exp(fl-s) 
{ ' k\ ~ {x-R)R k 

for any < R < x. 

3. Truncation of the L-series. 

To compute to D digits of accuracy, we need enough terms N in 
the L-series so that 



(7) 2^r^(n)|Re(G( S ,27rn/ g 1 / 2 ))| < 10' 



n>N 

Since the r^(n) are the coefficients of a weight one cusp form, r^(n) <C 
n 1/2 . Using ©, the left side of is 

< ^(g/n) 1/2 exp(-27m/g 1/2 ) 



n>N 

"OO 



/>oo 

~9 1/2 / y' 1/2 exp(-2ny/q 1/2 )dy 
Jn 

POO 

<{q/Nf' 2 / exp(-2ny/q^)dy 



JN 

<g/iV 1/2 exp(-27riV/g 1/2 ). 
For this to be less than 10~ D we need 

log(g) + Dlog(10) < 2nN/q 1/2 + log(JV)/2; 
it suffices that 

(8) iV = g 1/2 log(g- 10 D )/27r. 



jeffrey stopple 

4. Arrangement of the Taylor expansions. 
We seek to find intervals of radius Aj centered at Xj and bounding 



circles of radius Rj so that 



< Aj < Rj < xj. 



The tail of the Taylor expansion 

oo 

Y,G hk {s)(^n/q l / 2 -^-\ k 

k=B 

is bounded, via © by 



k=B x j ~ Rj x j ~ i — A j/Rj 

where Gj t k(s) is as in ||TJ|. 

One might wish to minimize the total number of terms in all the 
Taylor expansions, but this seems intractable. Instead we will simply 
choose the truncation parameter B to be independent of j. To make 
this happen we require that 

Aj/Rj— (constant) l/c<l, 
Xj — Rj = (constant) K > 0. 

Furthermore, in order that the endpoints of the intervals meet, we 
require that 

Xj + Aj = Xj+i — Aj_|_i. 

This gives 

Xj = Rj + K = cAj + K 
Xj+i = Rj+i + K = cA j+1 + K, 



so 



Thus 



Now 



Xj+1 - Xj = c(A j+1 - Aj) 
Xj+i Xj = Aj -f- Aj+i. 

A j+1 + A, = c(A j+1 - Aj) 
A J+1 =A,(c+l)/(c-l). 
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and 

implies that 



A? = Rj/c = (xj - K)/c 



x j+1 = Xj {c + l)/(c - 1) - 2K/(c - 1). 

For simplicity we choose the parameters c = 2 and = xi/2 which 
gives 

(10) x. 
and 

(11) A,- = 3 ; 



= y (3- 1 + 1) , 



The simplest choice for x 1 is 2n/q 1 ^ 2 . This gives 



4 = 



2tt 



1/2 



3 J 



-l 



2 3 J 



4 



Observe that I\ contains only the n = 1 term, J 2 contains only the 
n = 2 term, 7 3 contains only n = 3, . . . , 7, etc. Thus for the very 
small values of j where the intervals contain fewer than B terms 
one can do better computing each term instead of using the Taylor 
expansion. This would improve the efficiency about 40% for q near 
10 6 , 30% near 10 9 , and 20% near 10 20 . 

We can now determine how many intervals T are needed, since 
we need to compute terms in the series out to 

n = N = q 1 ' 2 log(g • 10 d )/2tt, 



so 



2nN 



\og(q ■ 10 D ) 



This requires that 



Xi 



x T + A T =^ Y + lj^> \og(q ■ 10 u ). 

With xi = 2n/q 1/2 we want 

2q 1/2 \og(q ■ 10 D )/n < 3 T , 

or 

T > {log(2/vr) + logfe 1 ' 2 • log(g • 10 D ))} / log(3). 

With 

log(2/7r)<0, and l/log(3)<l 
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we can simply choose 

(12) T = \ogiq 1 ' 2 ■ \og(q ■ 10 D )) 

5. Truncation of the Taylor expansions. 
Suppose we make error 5 in computing 

G(s, Zitn/q 1 ' 2 ) + G{1 - s, 2 7 m/g 1 / 2 ). 

Using the bound n 1 ^ 2 for r^(n) and the rough estimate q 1 ^ 2 for N, the 
maximal error we make in the sum for Z(t, (j>) is bounded by 

gl/2 

g- 1 / 4 ^n 1 / 2 -<5«5-g 1 / 2 . 

n=l 

Alternatively we can consider the standard error, assuming the er- 
rors in the terms are independent with standard deviation e. Then 
the standard error in the sum is bounded by |2 1 

(*» \ 1/2 

g" 1 / 4 X> 1/2 " e ) 2 « e ^ 1/4 - 

We will assume the latter from now on. We want e ■ g 1 / 4 < 1CT D , or 

(13) e < <r 1/4 l(T D , 

which will determine how many terms B we need to take in each 
Taylor expansion. 

It follows from the choice made in §3 that 

-J- — - and Xj — Rj = X\l2. 
Thus the tail © of the Taylor series reduces to 

(14) e = 2ex P(~ a; i/ 2 ) 2 i-i? 

X\ 

Again, with X\ = 2n /q 1 ^ 2 we combine (fT3l and <fl4|> to find that we 
require 

g 1 / 2 exp(-7r/g 1 / 2 )2 1 - B /7r < q-^lO' 13 . 

Since 

exp(-7r/g 1/2 )2/7r < 1 
we ignore it, and so we need g 3 / 4 10 D < 2 B , and let 

(15) B = 1.5 • log(g 3/4 10 D ) > log 2 (g 3/4 10 D ). 
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Even under the assumption of maximal error we would still only 
need 5 = O(log(gl0 D )). 

6. Implementation and examples. 

To implement this algorithm requires the computation of reduced 
representatives of all the forms, as well as their coordinates in terms 
of the generators of the cyclic factors of the class group. This com- 
putation is clearly 0(q 1 ^ 2 ). 

We then need to evaluate all the forms ax 2 + bxy + cy 2 on a rect- 
angular grid of integer lattice points (x, y), large enough not to miss 
any representation of any integer n < N. Determining the dimen- 
sions of the grid is a Lagrange multipliers problem; we must maxi- 
mize g(x,y) = x (respectively f(x,y) = y) subject to the constraint 
ax 2 + bxy + cy 2 = N. One finds 

x < 2(Nc/q) 1/2 and y < 2(Na/q) 1/2 . 

Since the forms are reduced, the inequalities on c and a and our 
choice 10! of N imply 

x < q 1/4 \og(q ■ 10 D ) 1/2 and y < (2/3 1/4 ) log(g ■ 10 D ) 1/2 . 

We have triple (x 2 , xy,y 2 ) for each lattice point, which are the rows 
of a 0(g 1//4 log(g • 10 D )) x 3 matrix. Evaluation of the all the forms at 
all the lattice points consists of multiplying the above matrix by the 
3 x 0{q 1 ^ 2 ) matrix of form data. This is 0(q 3 ^ 4 log(g)) multiplications 
(in terms of q, for a fixed number D of digits of accuracy.) 

The representation numbers r<g(n) are computed by brute force 
and ignorance; we look at each entry in the matrix product and in- 
crement the corresponding rg(ri), another 0(q 3 ^ log(g)) operations. 

Fast Fourier Transform on a group of size 0(q 1 ^ 2 ) is 0(q 1 ^ 2 log(g 1 / 2 )), 
and the function is vector valued with N = 0(q 1 ^ 2 log(g)) entries, for 
a total computation of size 0(q\og(q) 2 ). 

For each character 0, the precomputation of the sums 

^ ^Trn/q^-Xj)* 

takes N-B operations, which is 0(q 1 ^ 2 log(g) 2 ) by © and (fT5b . Subse- 
quently, evaluations of Z(t, 4>) using (J2J) require only T ■ B operations, 
which is Ottoeia) 2 ) by Ol and (Il5b. 

This algorithm was implemented in Mathematica and run on a 400 
MHz. Apple Powerbook G4 under OS X 10.2 1 . Mathematica is a good 

^This paper was written in June 2003. 
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TABLE 1. Class group data for various discriminants 



choice if one wants to re-invent the wheel as few times as possible. 
The incomplete Gamma function T(s, x) is available; it is computed 
via hypergeometric functions and continued fractions according to 
|9j A.9.4]. Also the Fast Fourier Transform is supported in a sophis- 
ticated way |9j A.9.4], via decomposition of the length of lists of data 
into prime factors. For large factors, fast convolution methods are 
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FIGURE 1. Histogram data of first zero above 0. 

used. Since our data are real (the representation numbers rg(n)), 
the algorithm makes use of a real transform method. Non-cyclic 
abelian groups are handled automatically. The Mathematica func- 
tion FindRoot finds zeros of functions by a combination of damped 
Newton's method, the secant method, and Brent's method |9j A.9.4]. 

On the other hand, a specialized package for number theory, such 
as PARI 1131 is preferable for computations involving binary qua- 
dratic forms. Fortunately, Mathematica 's MathLink capability allows 
users to install their own C code, and this was done, essentially in- 
stalling all the PARI functions for computations with binary qua- 
dratic forms into Mathematica. (This was documented in |14J.) 

The accuracy of the algorithm was checked by computing, to 15 
digits, the zeros of some genus character L-functions. By genus the- 
ory, these are products of Dirichlet L-functions attached to quadratic 
characters. These zeros were first computed directly, writing the 
Dirichlet L-function as a linear combination of Hurwitz zeta func- 
tions (also supported in Mathematica). The zeros were also compared 
to the data in |16[, and agreed to the number of digits given. 

The algorithm was then run on 21336 L-functions with conductor 
near q = 10 7 . The Hardy function Z(t, 4>) was evaluated, to six digits 
of accuracy, for t between 0. and 1. in steps of size 27r/(20 log(g)), 
looking for sign changes. When detected, Mathematica 's FindRoot 
was used to find the zero. This calculation took 137 hours, and found 
35190 zeros below t = 1. No attempt was made to prove GRH for 
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FIGURE 2. 1-level density vs. 1 - sin(27rx)/(27rx) 

these L-functions in this range, or to prove that all the zeros had been 
located. The zeros themselves are available at 
http : / / www . math . ucsb . edu/ ~stopple/ 

Table H] shows the discriminants, class numbers, class group struc- 
ture, and number of L-functions (omitting real (genus) characters, 
and one of each complex conjugate pair.) Figure [T] shows a histogram 
of of the lowest zero, renormalized by 7 = 7 • log(q) / (2ir) . These 
seem to support a conjecture that the corresponding random matrix 
model is symplectic (see Figure 5]). The mean height above 0. is 
1.13, considerably larger than 



0.78 



tvi(U Sp)(t)dt 



predicted by this model, due to the slow convergence as q 
Figure 12 compares the 1-level density of U Sp(oo), that is, 



00. 



1 — sin(27nr)/(27ra;) 



to the histogram data 



j^^££W0.78/1.18.7), 



where 



[a, (3) = [0, .05), [.05, 0.1), . . . , [1.75, 1.80) 
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FIGURE 3. Mean of first zero vs. h(-q)/^/q. 



The zeros 7 have been re-renormalized so that the mean above 0. of 
the first zero is 0.78. 

Finally, for each of the 29 discriminants considered, the mean above 
0. of the first zero for the L-i unctions with that discriminant was 
compared to the relative size h(—q) / y/q of the class number. Figure |3] 
shows the plot. The discriminants with smaller class number tend to 
have a larger mean first zero. The correlation between the two quan- 
tities is —0.89. The connection between this phenomenon and the 
Chowla-Selberg formula is discussed in §3 of [1J. 

As mentioned in the introduction, this method can be adapted to 
compute other types of L-functions. For an elliptic curve L-function, 
one gives up the advantage of the Fast Fourier Transform, but gains 
the advantage of doing the precomputations with integer arithmetic. 
The elliptic curve 

E : x 3 + y 3 = 763002 
has rank 5 according to |3|, so by the Birch-Swinnerton Dyer conjec- 
tures the L functions L(s, E) should have a zero of order 5 at s = 1. 
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FIGURE 4. Elliptic curve L-function with conjectural 
order 5 zero at s = 1 (i.e. t = 0). 

The conductor is 1746516156012. Figure H] show a plot of the corre- 
sponding Hardy function Z(t, E). 

Low lying zeros of L-functions for quadratic Dirichlet characters 
Xd are useful and scarce. The algorithm described here is of no real 
use since the existence of such a zero can be ruled out by a single 
evaluation. However, highly accurate values of zeros have been use- 
ful in the past, for example 1111 . With the parameter D = 110, and 
discriminant d = —175990483 we can compute L(l/2 + it, Xd) is zero, 

t = 0.0004752439954201629008767557526752 

446841851348862432434240449427326484 

62812721184470556544512670480839630. 

This was checked with the method of [17J, and is accurate to more 
than 100 digits. 
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